/***************************************************************
                       score1.c

only routines necessary for serial 
**************************************************************/
#include <stdio.h>
#include <malloc.h>
#include <float.h>
#include <math.h>
#include "clam.h"

#define FAST_SUL_ITERATIONS 3
#define graphQuery messQuery

void CpMdisplay(), ZautoCpM();
int  ZscoreCpM();
BOOL ZboolCpM();

static CLONE my_clone;
extern void scanPzA(), scanPzB(), scanPz2(), EvalCtgDisplay();
extern void CBCtgDisplay(), ClamCtgDisplay2();
extern void ZclearHigh(), ZhighCin();
extern int Zget_tol();
extern void ClamShare();
extern int NoPop;
extern int validDistribution;
extern float *distrib;

BOOL ZboolCpM();
int Zcnt_shared_markers();

/*************************************************************
                     CONTIGC
Finds the probability of coincidence
C(n, k) = (N over K) = bionomial coeffient = N!/K!(N-K)!

sum of N to nbandL C(nbandL, N) * (1 - p)^N * p^(nbandL-N)
where p = (1 - (2*tol/gel-len)) ^ nbandH

**************************************************************/
#define abs(x) ((x) >= 0 ? (x) : -(x))

/***********************************************************
                 DEF: Zmott
***********************************************************/
void Zmott()
{
double bin, b, lp1, lp2, p1, p2, p, tol, tol1, pt;

    p1 = (double) C1z.nbands / (double) Proj.gel_len;
    p2 = (double) C2z.nbands / (double) Proj.gel_len;
    p = p1*p2;
    lp1 = log(1.0-p1);
    lp2 = log(1.0-p2);
    tol = (double) Pz.tol;
    tol1 = (double) Pz.tol+1;
    pt = -p + p1*(1.0 - exp(tol1*lp2))*(exp(tol*lp2)) +
              p2*(1.0 - exp(tol1*lp1))*(exp(tol*lp1));
    if (validDistribution) {
      b = (double) Sz.newMatch;
    } else {
      b = (double) Sz.match;
    }

    bin = (double) Proj.gel_len;
    Sz.prob = (double) 
      exp(b*log(pt) + (bin-b)*log(1.0-pt) + 
      lgamma(bin+1.0) -lgamma(b+1.0) - lgamma(bin-b+1.0));
}
/***********************************************************
                 DEF: fast_sulston
*************************************************************/
static void fast_sulston()
{
register int i, k, n, nbH, nbL;
double a, c, pp, pa, pb;
static int first_time=1, tol = -1;
static double lim, Logfact[NBANDS+1], psmn=0.0;

    if (first_time) {
       Logfact[0] = 0.0;
       Logfact[1] = 0.0;
       for (i = 2; i <= NBANDS; ++i)
         Logfact[i] = Logfact[i - 1] + log((double) i);
       first_time=0;
       lim = log(DBL_MIN);
    }
    if (Pz.tol == 0)  psmn = 1 - (double) 1.0 / Proj.gel_len;
    else  psmn = 1 - (double) (Pz.tol << 1) / Proj.gel_len;
    tol = Pz.tol;

    nbL = MiN(C1z.nbands,C2z.nbands);
    nbH = MaX(C1z.nbands,C2z.nbands);

    Sz.prob = DBL_MIN;
    
    pp = (double) pow(psmn, nbH);
    pa = log(pp);
    pb = log(1.0 - pp);

    /* 18june00 starting at nbl 64, nhl 74, if the match is zero,
            this optimization causes the prob to be below cutoff */

    if (Sz.match==0) Sz.prob =1.0;
    else {
      for (k=0, n = Sz.match; n <= nbL && k < FAST_SUL_ITERATIONS; n++, k++)
      {
            c = Logfact[nbL] - Logfact[nbL - n] - Logfact[n];
            a = c + n * pb + (nbL - n) * pa;
            Sz.prob +=  (a >= lim) ? exp(a) : exp(lim);
       }
    }
}
/***********************************************************
                 DEF: pure_sulston
*************************************************************/
static void pure_sulston()
{
register int i, n, nbH, nbL;
double a, c, pp, pa, pb;
static int first_time=1;
static double lim,  Logfact[NBANDS+1], psmn=0.0;

    if (first_time) {
       Logfact[0] = 0.0;
       Logfact[1] = 0.0;
       for (i = 2; i <= NBANDS; ++i)
         Logfact[i] = Logfact[i - 1] + log((double) i);
       first_time=0;
       lim = log(DBL_MIN);
    }
    nbL = MiN(C1z.nbands,C2z.nbands);
    nbH = MaX(C1z.nbands,C2z.nbands);

    Sz.prob = DBL_MIN;
    
    pp = (double) pow(psmn, nbH);
    pa = log(pp);
    pb = log(1.0 - pp);

    for (n = Sz.match; n <= nbL; n++)
    {
       c = Logfact[nbL] - Logfact[nbL - n] - Logfact[n];
       a = c + n * pb + (nbL - n) * pa;
       Sz.prob +=  (a >= lim) ? exp(a) : exp(lim);
    }
}
/***********************************************************
                 DEF: Zsulston
*************************************************************/
void Zsulston(int min)
{
    Zgetmatch();
    if (Proj.eq2==1) Zmott();
    else if (Proj.eq2==2 || C1z.nbands >=60 || C2z.nbands >= 60) 
         pure_sulston();
    else 
         fast_sulston();
}

/*********************************************************************
                    DEF: Zgetmatch

This accepts two clones and computes the number of matching bands
which it stores in result.
*******************************************************************/
void Zgetmatch()
{
  int kbnd, k, l, idiff, lstart;

  Sz.total = Sz.match = Sz.newMatch = 0;
  lstart = 0;
  for (k = 0; k < C1z.nbands; ++k)
  {
      kbnd = C1z.coords[k];
      for (l = lstart; l < C2z.nbands; ++l)
      {
          idiff = kbnd - C2z.coords[l];
          if (Ztol(kbnd, idiff))
            {
              Sz.match++;
              if (validDistribution) {
                /* trying different values for distribution */
                Sz.newMatch += (distrib[kbnd] + distrib[kbnd-idiff])/2;
                Sz.match = Sz.newMatch;
              }
              lstart = l + 1;
              break;
            }
          else if (idiff < 0)
            {
              lstart = l;
              break;
            }
       }
   }
  return;
}
/**************************************************************
            DE: ZscoreCpM
All clam routines call this for scoreing (Zsulston or Zmott must be called
before this).
**************************************************************/
int ZscoreCpM()
{
int junk, exponent=0;
char str[20];
register int i;

   Sz.mark = 0;
   if (Sz.prob < Pz.cutFlt) {
       CpMTbl[0].cnt++;
       sprintf(str,"%.0e", Sz.prob);
       sscanf(str,"%de-%d",&junk, &exponent);
       if (exponent==0) {
           printf("Clone1 %i\tScore%.0e\tScore%f\n", 
                     C1z.nbands,Sz.prob,Sz.prob);
           printf("exp = 0\n");
       }
       return exponent;
   }
   if (Cp.useFlag==0)  return 0;
   if (!Zcnt_shared_markers()) return 0;
   if (Sz.prob > CpM[MAX_CpM-1].cut) return 0;

   for (i=1; i<MAX_CpM; i++)
      if (Sz.prob < CpM[i].cut && Sz.mark >= CpM[i].nm) {
         CpMTbl[i].cnt++;
         sprintf(str,"%.0e", Sz.prob);
         sscanf(str,"%de-%d",&junk, &exponent);
         if (exponent==0) printf("exp = 0\n");
         return exponent;
      }
   return 0;
}
/******************************************************************
                 DEF: Zcnt_shared_markers()
*****************************************************************/
int Zcnt_shared_markers()
{
struct marker *markerptr;
struct markertop *mk1, *mk2;
CLONE *c1, *c2;
int score=0;

   c1 = arrp(acedata, C1z.cin, CLONE);
   c2 = arrp(acedata, C2z.cin, CLONE);
   for (mk1 = c1->marker; mk1 != NULL; mk1 = mk1->nextmarker)
   {
        markerptr = arrp(markerdata, mk1->markerindex, MARKER);
        if(!Cp.useYBP &&
             (markerptr->type == markYAC || markerptr->type == markBAC ||
                                markerptr->type == markPAC)) continue;
        if(!Cp.usePCR && markerptr->type == markPCR) continue;
        if(!Cp.useREP && markerptr->type == markREP) continue;

        for (mk2 = c2->marker; mk2 != NULL; mk2 = mk2->nextmarker)
            if (mk1->markerindex == mk2->markerindex) {
               score++;
               break;
            }
   }
   Sz.mark = score;
return score;
}

